library(readxl) library(readODS) library(writexl) library(ggthemes) library(lvplot) library(glue) library(tidyverse) library(plotly) # gráfico 3d library(randtests) library(goft)
Tópico 5
Pacotes que iremos usar na semana 6
Gerando número aleatórios
- Geralmente, números aleatórios são gerados por algoritmos matemáticos
- Algoritmos matemáticos são determinísticos, e usamos o termo pseudoaleatório
- Amostra construída com números pseudoaleatórios: pseudo-amostra
- Geração de números aleatórios de distribuições de probabilidade é baseada na geração de números aleatórios da distribuição uniform
- É possível coletar valores aleatórios da distribuição uniforme que verdadeiramente aleatório usando fenômenos naturais (como ruídos atmosféricos e o magnetismo da terra).
- Em geral, evitamos o uso de números verdadeiramente aleatórios em favor dos número pseudoaleatórios: números pseudoaleatórios são reprodutíveis
- Números pseudoaleatórios são baseados em um valor chamado chave ou semente
- Uma sequência de números pseudoaleatórios geralmente é periódica
Gerador de números pseudoaleatórios
Definição 1 (Gerador de números pseudoaleatórios) Um número gerador pseudoaleatório é uma estrutura \(\Xi = (S, s_0, T, U, G)\), \(S\) é um conjunto de números (chamados estados), \(s_0\) é a chave ou semente, \(T:S \rightarrow S\) é uma função (chamada de função de transformação), \(U\) é um conjunto de números (chamados símbolos de saída), e \(G:S \rightarrow U\) é uma função (chamada de função de saída).
- \(s_n = T(s_{n-1}), n=1, 2, 3, \dots\)
- Como \(\lvert S \rvert < \infty\), então existe \(0 \leq i < j\) tal que \(s_i = s_j\) e \(s_{i+n} = s_{j+n}, \forall n \geq 0\)
- \(s_{i + n} = s_{j + n} \Rightarrow s_{i + n} = s_{i + (j - i) + n}, \forall n \geq 0 \Rightarrow s_{m} \geq s_{(j - i) +m}, \forall m \geq i\)
- Período: menor inteiro \(p\) tal que \(s_{n + p} = s_{n}\) para \(n \geq r\) para algum \(r \geq 0\)
- Note que \(p \leq \lvert S \rvert\)
- Desejável: \(p \approx \lvert S \rvert\) (poucas repetições na sequência de números pseudoaleatórios)
Gerador congruente linear
Definição 2 (Gerador Congruente Linear) Seja \(m>0\), \(0<a<m\), \(0 \leq c \leq m\) e \(0 \leq s_0 < m\). Considere o gerador pseudoaleatório \(\Xi = (S, s_0, T, U, G)\) onde \[T(s_i) = (a \cdot s_{i-1} + c) \mod m, \quad 0 \leq s_i < m, \quad i=0, 1, 2, 3, \dots.\] Chamamos \(a\) de multiplicador, \(m\) de módulo, \(c\) de incremento, \(s_0\) de chave ou semente e \(\Xi\) é denominado de Gerador Congruente Linear.
O gerador congruente linear mais famoso é o randu proposto em 1951 e implementado pela IBM logo em seguida.
RANDU
Definição 3 (RANDU) Considere o Gerador Congruente Linear \(\Xi\), em que \(S = \mathbb{Z} \cap [0, 2^{31}]\), \(s_0 \in S\), \(T(s) = (2^{16} + 3) s \mod 2^{31}\), \(U = \frac{S}{2^{31}}\) e \(G(s) = \frac{s}{2^{31}}\).
Restrições
Suponha que \(s_{k+1} < 2^{31}\), então \[ \begin{split} s_{k+2} &= (2^{16} + 3) s_{k+1} = (2^{16} + 3) (2^{16} + 3) s_{k},\\ &= (2^{32} + 6 \cdot 2^{16} + 9) s_k = 6 (2^{16} + 3) s_k - 9 s_k,\\ &= 6 s_{k+1} - 9 s_k. \end{split} \]
Ou seja, existe uma relaçãos simples entre \(x_{k+2}\), \(x_{k+1}\) e \(x_k\), e podemos descobrir o gerador pseudoaleatório a partir da sequência de número pseudoaleatórios.
RANDU
randu <- function(n, seed = 1) {
x <- vector(mode = "double", n)
a <- 2^16 + 3
m <- 2^31
x[1] <- (a * seed) %% m
for (i in 2:n) x[i] <- (a * x[i - 1]) %% m
x / m
}
Marsaglia (1968) provou os valores desta matriz estão dentro de 15 hiperplanos.
Vamos usar o teste de aleatoriedade (Wald e Wolfowitz 1940): runs.test() do pacote randtests.
Criando uma amostra com RANDU
tamanho_amostra <- 9999 amostra <- randu(tamanho_amostra) matriz <- tibble( xk = amostra[seq_len(tamanho_amostra) %% 3 == 0], xk1 = amostra[seq_len(tamanho_amostra) %% 3 == 1], xk2 = amostra[seq_len(tamanho_amostra) %% 3 == 2] )
Checando a amostra
ks.test(amostra, "punif")
## ## One-sample Kolmogorov-Smirnov test ## ## data: amostra ## D = 0.0063368, p-value = 0.8168 ## alternative hypothesis: two-sided
runs.test(amostra)
## ## Runs Test ## ## data: amostra ## statistic = -2.0203, runs = 4899, n1 = 4999, n2 = 4999, n = 9998, p-value = 0.04335 ## alternative hypothesis: nonrandomness
Gráfico com RANDU
plot_ly(matriz, x = ~xk, y = ~xk1, z = ~xk2, size = 0.2)
Gerador Defasado de Fibonacci
Definição 4 (Gerador Defasado de Fibonacci) Seja \(i,j,k\) números inteiros tal que \(0 < k < j < i\), e considere o Gerador Congruente Linear \(\Xi\) com em que \(S = \mathbb{Z} \cap [0, 2^{31}]\), \(T(s_i) = (s_{i-j} + s_{j-k}) \mod 2^{31}\), \(U = \frac{S}{2^{31}}\), \(G(s) = \frac{s}{2^{31}}\) e \(s_0\) é um conjunto de valores iniciais tal que \(\lvert s_0 \rvert > j\). Chamamos \(\Xi\) de Gerador Desafado de Fibonacci.
Este algoritmo é uma melhoria sobre o randu, e pode ser usado para simulações. (Mas não é recomendado para criptografia).
O R implementou o algoritmo Mersenne Twister (vide Härdle, Okhrin, e Okhrin 2017 para mais detalhes) com a função runif(n).
GDF
k <- 5
j <- 17
m <- 2^31
semente <- runif(j, 0, m) |> round()
gdf <- function(n, seed = semente) {
x <- seed
for (i in seq_len(n)) {
x[j + i] <- (x[i] + x[j + i - k]) %% m
}
x[(j + 1):length(x)] / m
}
Amostra
tamanho_amostra <- 9999 amostra <- gdf(tamanho_amostra) matriz <- tibble( xk = amostra[seq_len(tamanho_amostra) %% 3 == 0], xk1 = amostra[seq_len(tamanho_amostra) %% 3 == 1], xk2 = amostra[seq_len(tamanho_amostra) %% 3 == 2] )
Checando a amostra
ks.test(amostra, "punif")
## ## One-sample Kolmogorov-Smirnov test ## ## data: amostra ## D = 0.0082316, p-value = 0.507 ## alternative hypothesis: two-sided
runs.test(amostra)
## ## Runs Test ## ## data: amostra ## statistic = 1.4002, runs = 5070, n1 = 4999, n2 = 4999, n = 9998, p-value = 0.1615 ## alternative hypothesis: nonrandomness
Gráfico com GDF
plot_ly(matriz, x = ~xk, y = ~xk1, z = ~xk2, size = 0.2)
Métodos para geração de outras distribuições
Método da Transformação Inversa
Definição 5 (Inversa Generalizada) Seja \(F\) uma função de distribuição acumulada, chamamos \(F^{-}(u) = \inf \{x \mid F(x) \geq u\}, \forall u \in [0, 1]\) de função inversa generalizada.
Proposição 1 (Método da transformação inversa.) Se \(U \sim U(0, 1)\), então \(F^{-}(U)\) é uma variável aleatória com FDA \(F\).
Método da Transformação Inversa
Prova. Note que
- \(F(z) \geq u \Rightarrow z \in \{x \mid F(x) \geq u\} \Rightarrow z \geq F^{-}(u) \Rightarrow [F(z) \geq U] \subset [z \geq F^{-}(U)]\)
- \(F^{-}(u) \leq z \Rightarrow u \leq F(z) \Rightarrow [F^{-}(U) \leq z] \subset [U \leq F(z)]\)
Método da Transformação Inversa
Algoritmo
- Encontre \(F^{-}(u)\)
- Gere um valor \(u\) de \(U \sim U(0, 1)\)
- Calcule \(F^{-}(u)\)
Como encontrar \(F^{-}\)
- Se \(X\) é uma variável aleatória contínua: \(F^{-1} = F^{-}\)
- Se \(X\) é uma variável alatória discreta com suporte \(\chi = \{x_1, \dots, x_k, \dots\}\): \(F^{-}(u) = \begin{cases} x_i, & F(x_{i-1}) < u \leq F(x_i), \\ \min(\chi), & u \leq F(\min(\chi)) \end{cases}\)
Método da Transformação Inversa
Exemplo (Caso Discreto):
- Seja \(X\) uma variável aleatória discreta com suporte \(\chi=\{0, 1, 2\}\) com função de probabilidade \(f(0)=0,3\), \(f(1)=0,2\) e \(f(2) = 0,5\).
- \(F(x) = \begin{cases} 0, & x < 0,\\ 0,3, & 0 \leq x < 1,\\ 0,5, & 1 \leq x < 2,\\ 1, & x \geq 2 \end{cases}\)
- \(F^{-}(u) = \begin{cases} 0, & u \leq 0,3,\\ 1, & 0,3 < u \leq 0,5\\ 2, & 0,5 < u \leq 1 \end{cases}\)
Método da Transformação Inversa
Exemplo (Caso Discreto):
ex_discreto <- function(n) {
if (length(n) > 1) stop("n precisa ser um valor escalar inteiro.")
seq_len(n) |>
map_dbl(\(x) {
u <- runif(1)
if(u <= 0.3) {
return(0)
} else if (0.3 < u && u <= 0.5) {
return(1)
} else if (u > 0.5) {
return(2)
}
})
}
Método da Transformação Inversa
Exemplo (Caso Discreto):
set.seed(121343) tibble(amostra = ex_discreto(1000)) |> group_by(amostra) |> summarise(freq = n()) |> mutate(fr = freq / sum(freq))
## # A tibble: 3 × 3 ## amostra freq fr ## <dbl> <int> <dbl> ## 1 0 304 0.304 ## 2 1 195 0.195 ## 3 2 501 0.501
Método da Transformação Inversa
Exemplo (Caso Contínuo):
- Seja \(X \sim Exp(\lambda)\) com função densidade \(f(x) = \lambda \exp(-\lambda x), x> 0\) e \(F(x) = 1 - \exp(-\lambda x), x> 0\)
- \(F^{-}(u) = \begin{cases} \frac{-\ln(1 - u)}{\lambda}, & u > 0,\\ 0, & u \leq 0, \end{cases}\)
lambda <- 2
exp_2 <- function(n) {
if (length(n) > 1) stop("n precisa ser um valor escalar inteiro.")
seq_len(n) |>
map_dbl(\(x) {
u <- runif(1)
if (u <= 0) {
return(0)
} else {
return(-log(1 - u) / lambda)
}
})
}
Método da Transformação Inversa
Exemplo (Caso Contínuo):
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = exp_2(n)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "Amostra", y = "Densidade de Frequência")
Método da Transformação Inversa
Exemplo (Caso Contínuo):
ks.test(dados$amostra, pexp, rate = 2)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.03783, p-value = 0.1143 ## alternative hypothesis: two-sided
Método de Aceitação-Rejeição
Ideia
- Difícil gerar valores de \(X \sim f\)
- Simples gerar valores de \(Y \sim g\)
- \(X\) e \(Y\) têm o mesmo suporte
- \(g(x)\) cobre \(f(x)\): \(\exists c>0\) tal que \(f(x) \leq c \cdot g(x)\)
- Gerar valores de \(X\) a partir de \(Y\)
Método de Aceitação-Rejeição
Algoritmo
- Encontre \(Y\) e \(c>0\) tal que \(f(x) \leq c \cdot g(x)\)
- Gere um valor \(y\) de \(Y\)
- Gere um valor \(u\) de \(U \sim U(0, 1)\)
- Se \(u \leq \frac{f(y)}{c\cdot g(y)}\), então aceitamos \(y\) como valor gerado de \(X\). Se \(u > \frac{f(y)}{c\cdot g(y)}\), repetimos 2. e 3.
Geralmente, usamos \(c\) como sendo o valor máximo de \(\frac{f(x)}{g(x)}\).
Método de Aceitação-Rejeição
Exemplo (Caso Contínuo)
- \(X \sim B(2, 1)\), \(\chi = (0,1)\), \(f(x) = 2 \cdot x, x \in (0, 1)\) e \(F(x) = x^2, x \in (0, 1)\)
- Vamos considerar \(Y \sim U(0, 1)\) com \(g(y) = x, x \in (0, 1)\)
- \(\frac{f(x)}{g(x)} = 2 = c\)
beta_21 <- function(n) {
if (length(n) > 1) stop("'n' precisa ser um escalar inteiro.")
seq_len(n) |>
map_dbl(\(i) {
y <- runif(1); u <- runif(1)
while(u > dbeta(y, shape1 = 2, shape2 = 1) / (2 * dunif(y))) {
y <- runif(1); u <- runif(1)
}
return(y)
})
}
Método de Aceitação-Rejeição
Exemplo (Caso Contínuo)
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = beta_21(n)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "Beta(2, 1)", y = "Densidade de Frequência")
Método de Aceitação-Rejeição
Exemplo (Caso Contínuo)
ks.test(dados$amostra, pbeta, shape1 = 2, shape2 = 1)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.030908, p-value = 0.295 ## alternative hypothesis: two-sided
Método de Aceitação-Rejeição
Exemplo (Caso Discreto)
- Seja \(X\) com \(\chi=\{1, 2, 3, 4, 5\}\) com
| \(x\) | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|
| \(f(x)\) | \(0,15\) | \(0,22\) | \(0,33\) | \(0,10\) | \(0,20\) |
- Seja \(Y\) com \(\chi = \{1, 2, 3, 4, 5\}\) com \(g(x) = 0,20, \quad x \in \chi\)
- O valor máximo do conjunto \(\{\frac{f(x)}{g(x)}\mid x \in \chi \}\) é \(c = 1,65\)
Método de Aceitação-Rejeição
Algoritmo – caso discreto
fp <- c(0.15, 0.22, 0.33, 0.10, 0.20)
const <- 1.65
mar_discreto <- function(n) {
if (length(n) > 1) {
stop("'n' precisa um escalar inteiro.")
}
seq_len(n) |>
map_dbl(\(i) {
y <- sample.int(5, 1); u <- runif(1)
while(u > fp[y] / (0.2 * const)) {
y <- sample.int(5, 1); u <- runif(1)
}
return(y)
})
}
Método de Aceitação-Rejeição
Algoritmo – caso discreto
n <- 1000 dados <- tibble(amostra = mar_discreto(n)) print(fp)
## [1] 0.15 0.22 0.33 0.10 0.20
dados |> group_by(amostra) |> summarise(freq = n()) |> mutate(fr = freq / sum(freq))
## # A tibble: 5 × 3 ## amostra freq fr ## <dbl> <int> <dbl> ## 1 1 146 0.146 ## 2 2 211 0.211 ## 3 3 340 0.34 ## 4 4 98 0.098 ## 5 5 205 0.205
Amostrando das principais distribuições de probabilidade
Distribuição normal
Algoritmo Box-Müller
Sejam \(U_1\sim U(0,1)\) e \(U_2\sim U(0,1)\) duas variáveis aleatórias contínuas independentes, e considere \[ \begin{split} X_1 &= \sqrt{-2\ln(U_1)}\cos(2\pi U_2),\\ X_2 &= \sqrt{-2\ln(U_1)}\sin(2\pi U_2), \end{split} \] então \(X_1\) e \(X_2\) são independentes, e \(X_1 \sim N(0,1)\) e \(X_2 \sim N(0,1)\).
- Neave (1973) provou que o algoritmo de Box-Müller quando usado com randu não tem distribuição normal padrão.
Distribuição normal
Algoritmo Box-Müller
box_muller <- function(n) {
if (length(n) > 1) stop("'n' precisa ser um escalar inteiro.")
seq_len(n) |>
map_dbl(\(i) {
u1 <- runif(1); u2 <- runif(1)
sqrt(-2 * log(u1)) * cos(2 * pi * u2)
})
}
Distribuição normal
Exemplo – Box-Müller
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = box_muller(n)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "Distribuição normal.", y = "Densidade de Frequência")
Distribuição normal
Exemplo – Box-Müller
ks.test(dados$amostra, pnorm)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.022602, p-value = 0.6866 ## alternative hypothesis: two-sided
Distribuição normal
Método polar
Sejam \(U_1\sim U(-1,1)\) e \(U_2\sim U(-1,1)\) duas variáveis aleatórias contínuas independentes. Suponha que \(\omega = U_1^2 + U_2^2 \leq 1\), e considere \[ \begin{split} X_1 &= \sqrt{\frac{-2\ln(\omega)}{\omega}} U_1,\\ X_2 &= \sqrt{\frac{-2\ln(\omega)}{\omega}} U_2. \end{split} \] Então é possível provar que \(X_1\) e \(X_2\) são duas variáveis aleatórias independentes, e \(X_1 \sim N(0, 1)\) e \(X_2 \sim N(0,1)\).
Para gerar valores de \(U(-1, 1)\) usamos o método da transformação inversa. Note que \(2 \cdot U -1 \sim U(-1, 1)\) com \(U \sim U(0,1)\).
Distribuição normal
Exemplo – método polar
metodo_polar <- function(n) {
if (length(n) > 1) stop("'n' precisa ser um escalar inteiro.")
seq_len(n) |>
map_dbl(\(i) {
u1 <- 2 * runif(1) - 1; u2 <- 2 * runif(1) - 1; w <- u1^2 + u2^2
while(w > 1) {
u1 <- 2 * runif(1) - 1; u2 <- 2 * runif(1) - 1; w <- u1^2 + u2^2
}
return(sqrt(-2 * log(w) / w) * u1)
})
}
Distribuição normal
Exemplo – método polar
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = metodo_polar(n)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "Valores amostrados", y = "Densidade de Frequência")
Distribuição normal
Exemplo – método polar
ks.test(dados$amostra, pnorm)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.021913, p-value = 0.7229 ## alternative hypothesis: two-sided
Distribuição gamma
- \(X \sim \Gamma(\alpha, \beta)\) com função densidade \(f(x) = \frac{1}{\beta\cdot \Gamma(\alpha)} x^{\alpha-1}\exp\left(-\frac{x}{\beta}\right)\)
- Se \(Y \sim \Gamma(\alpha, 1)\), então \(\beta\cdot Y \sim \Gamma(\alpha, \beta)\)
Ideia
- Algoritmos para \(0 < \alpha \leq 1\)
- Algoritmos para \(\alpha > 1\)
Distribuição gamma
\(\alpha > 1\): proposto por Cheng (1977)
- Passo 0) Calcule \(a = \sqrt{\frac{1}{2\alpha -1}}\), \(b = \alpha - \ln(4)\) e \(c = \alpha + \frac{1}{a}\)
- Passo 1) Gere valores \(u_1 \sim U(0,1)\) e \(u_2 \sim U(0,1)\)
- Passo 2) Calcule \(v = a \cdot \log\left( \frac{u_1}{1-u_1} \right)\) e \(x = \alpha \exp(v)\)
- Passo 3) Se \(b + c v - x \geq \ln(u_1^2 u_2)\), então aceite que \(x\) é um valor gerado de \(\Gamma(\alpha, 1)\). Caso contrário, repita passos 1) e 2)
- Passo 4) Calcule \(\beta \cdot x\)
Distribuição gamma
\(0 < \alpha \leq 1\): proposto por Ahrens e Dieter (1974)
- Passo 0) Calcule \(b = \frac{\alpha + e}{e}\), em que \(e=\exp(1)\)
- Passo 1) Gere um valor \(u_1 \sim U(0,1)\) e calcule \(p = b \cdot u_1\)
- Passo 2) Gere um valor \(u_2 \sim U(0,1)\) e calcule \(E = -\ln(u_2)\)
- Passo 3)
- Se \(p \leq 1\), então calcule \(x = p^{\frac{1}{\alpha}}\). Se \(x \leq E\), então aceite que \(x\) é um valor gerado de \(\Gamma(\alpha, 1)\). Caso contrário, repita os passo 1) e 2).
- Se \(p > 1\), então calcule \(x = - \ln\left( \frac{b - p}{\alpha} \right)\). Se \((1-\alpha)\ln(x) \leq x\), então aceite que \(x\) é um valor gerado de \(\Gamma(\alpha, 1)\). Caso contrário, repita os passos 1) e 2).
Distribuição gamma
Exemplo: \(\alpha > 1\)
gamma_alpha_grande <- function(n, shape, scale) {
if (length(n) > 1) stop("'n' precisa ser um escalar inteiro.")
if (shape <= 1) stop("'shape' precisa ser um escalar maior que 1.")
if (scale < 0) stop("'scale' precisa ser um escalar positivo.")
seq_len(n) |>
map_dbl(\(i) {
a <- 1 / sqrt(2 * shape - 1); b <- shape - log(4)
ce <- shape + 1 / a
u1 <- runif(1); u2 <- runif(1)
v <- a * log(u1 / (1 - u1)); x <- shape * exp(v)
while(b + ce * v - x < log(u1^2 * u2)) {
u1 <- runif(1); u2 <- runif(1)
v <- a * log(u1 / (1 - u1)); x <- shape * exp(v)
}
return(scale * x)
})
}
Distribuição gamma
Exemplo: \(\alpha > 1\)
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = gamma_alpha_grande(n, shape = 2, scale = 4)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "valores amostrados", y = "Densidade de Frequência")
Distribuição gamma
Exemplo: \(\alpha > 1\)
ks.test(dados$amostra, pgamma, shape = 2, scale = 4)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.020675, p-value = 0.7861 ## alternative hypothesis: two-sided
Distribuição gamma
Exemplo: \(0 < \alpha \leq 1\)
gamma_alpha_pequeno <- function(n, shape, scale) {
seq_len(n) |>
map_dbl(\(i) {
b <- (shape + exp(1)) / exp(1); u1 <- runif(1); u2 <- runif(1); p <- b * u1; E <- -log(u2)
laco <- TRUE
while(laco) {
if (p <= 1) {
x <- p^(1 / shape)
if (x > E) {
b <- (shape + exp(1)) / exp(1); u1 <- runif(1); u2 <- runif(1); p <- b * u1; E <- -log(u2)
} else {
laco <- FALSE
}
} else {
x <- -log((b - p) / shape)
if ((1 - shape) * log(x) > E) {
b <- (shape + exp(1)) / exp(1); u1 <- runif(1); u2 <- runif(1); p <- b * u1; E <- -log(u2)
} else {
laco <- FALSE
}
}
}
return(scale * x)
})}
Distribuição gamma
Exemplo: \(0 < \alpha \leq 1\)
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = gamma_alpha_pequeno(1000, shape = 0.5, scale = 4)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "valores amostrados", y = "Densidade de Frequência")
Distribuição gamma
Exemplo: \(0 < \alpha \leq 1\)
ks.test(dados$amostra, pgamma, shape = 0.5, scale = 4)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.024642, p-value = 0.5782 ## alternative hypothesis: two-sided
Distribuição beta (Atkinson e Pearce 1976)
Seja \(X \sim B(\alpha, \beta)\) com função densidade dada por \(f(x) = \frac{\Gamma(\alpha) \Gamma(b)}{\Gamma(\alpha + \beta)} x^{\alpha - 1} (1 - x)^{\beta -1}\).
- Passo 1) Gere valores \(u_1 \sim U(0,1)\) e \(u_2 \sim U(0,1)\)
- Passo 2) Calcule \(v_1 = u_1^{\frac{1}{\alpha}}\), \(v_2 = u_2^{\frac{1}{\beta}}\) e \(w = v_1 + v_2\)
- Passo 3) Se \(w \leq 1\), então aceite que \(x\) é um valor gerado de \(B(\alpha, \beta)\). Caso contrário, repita os passos 1) e 2).
Distribuição beta (Atkinson e Pearce 1976)
Exemplo
r_beta <- function(n, shape1, shape2) {
seq_len(n) |>
map_dbl(\(i) {
u1 <- runif(1); u2 <- runif(1); v1 <- u1^(1 / shape1)
v2 <- u2^(1 / shape2); w <- v1 + v2
while(w > 1) {
u1 <- runif(1); u2 <- runif(1); v1 <- u1^(1 / shape1)
v2 <- u2^(1 / shape2); w <- v1 + v2
}
return(v1 / w)
})
}
Distribuição beta (Atkinson e Pearce 1976)
Exemplo
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = r_beta(n, shape1 = 2, shape2 = 3)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "valores amostrados", y = "Densidade de Frequência")
Distribuição beta (Atkinson e Pearce 1976)
Exemplo
ks.test(dados$amostra, pbeta, shape1 = 2, shape2 = 3)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.035974, p-value = 0.1502 ## alternative hypothesis: two-sided
Distribuição binomial
Seja \(X \sim b(n, p)\).
- Passo 0) Defina inicialmente \(x = 0\)
- Passo 1) Gere um valor \(u \sim U(0,1)\)
- Passo 2) Se \(u \leq p\), então acrescente \(1\) a x
- Passo 3) Repita os passos 1) e 2) \(n\) vezes
Distribuição binomial
Exemplo
r_binom <- function(n, size, p) {
seq_len(n) |>
map_dbl(\(i) {
seq_len(size) |> map_dbl(~ runif(1) <= p) |> sum()
})
}
Distribuição binomial
Exemplo
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = r_binom(1000, 5, 0.75)) ggplot(data = dados) + geom_bar(aes(x = amostra, y = 100 * ..prop..)) + theme_minimal() + labs(x = "valores amostrados", y = "Porcentagem")
Distribuição binomial
Exemplo
n <- 1000 size <- 5; p <- 0.55 suporte <- 0:size fda <- stepfun(x = suporte, y = c(0, pbinom(suporte, size = size, prob = p))) dgof::ks.test(rbinom(1000, size, p), fda)
## ## One-sample Kolmogorov-Smirnov test ## ## data: rbinom(1000, size, p) ## D = 0.024218, p-value = 0.6006 ## alternative hypothesis: two-sided
Distribuição poison
Seja \(X \sim Poison(\lambda)\).
- Passo 0) Defina inicialmente \(k = 0\) e \(p = 1\)
- Passo 1) Acrescente \(1\) a \(k\), gere um valor de \(u \sim U(0, 1)\) e atualize \(p\) para \(p\cdot u\)
- Passo 2) Se \(\exp(-\lambda) >= p\), aceite \(k-1\) como um valor gerado da distribuição \(Poison(\lambda)\). Caso contrário, repita 1).
Distribuição poison
Exemplo
r_pois <- function(n, lambda) {
seq_len(n) |>
map_dbl(\(i) {
k <- 0; p <- 1
while(exp(-lambda) < p) {
k <- k + 1; p <- p * runif(1)
}
return(k - 1)
})
}
Distribuição poison
Exemplo
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = r_pois(1000, lambda = 4)) ggplot(data = dados) + geom_bar(aes(x = amostra, y = 100 * ..prop..)) + theme_minimal() + labs(x = "valores amostrados", y = "Porcentagem")
Distribuição poison
Exemplo
suporte <- 0:10000 fda <- stepfun(suporte, c(0, ppois(suporte, lambda = 4))) dgof::ks.test(dados$amostra, fda)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.026163, p-value = 0.5004 ## alternative hypothesis: two-sided
Distribuição multivariada
Distribuição normal multivariada
- Seja \(\mathbf{X}\) um vetor aleatório. Lembre que \(E(A\cdot \mathbf{X} + \mathbf{b}) = A\cdot E(\mathbf{X}) + b\), e \(Var(A\cdot \mathbf{X} + \mathbf{b}) = A\cdot Var(\mathbf{X})\cdot A^\top\)
- \(Se \mathbf{X} \sim N(\mathbf{\mu}, \Sigma)\), então \(R^{-1}(\mathbf{X} - \mathbf{\mu}) \sim N(\mathbf{0}, I)\) em que \(R\) é decomposição de Cholesky (\(R\cdot R^\top = \Sigma\))
- \(Se \mathbf{X} \sim N(\mathbf{0}, \Sigma)\), então \(R\cdot \mathbf{X} + \mathbf{\mu} \sim N(\mu, \Sigma)\)
Exemplo
r_mnorm <- function(n, sigma, media) {
media <- matrix(media, ncol = 1)
R <- chol(m_var)
seq_len(n) |>
sapply(\(i){
return(R %*% matrix(nrow(sigma) |> rnorm(), ncol = 1) + media)
}) |>
t()
}
Distribuição normal multivariada
Exemplo
p <- 0.90
dim <- 10
m_var <- seq_len(dim) |>
sapply(\(i) {
seq_len(dim) |> map_dbl(\(j) p^abs(i - j))
})
v_media <- rep(2, dim)
dados <- r_mnorm(1000, sigma = m_var, media = v_media)
goft::mvshapiro_test(dados)
## ## Generalized Shapiro-Wilk test for Multivariate Normality ## ## data: dados ## MVW = 0.99832, p-value = 0.5471
Referências
Ahrens, Joachim H, e Ulrich Dieter. 1974. «Computer methods for sampling from gamma, beta, poisson and bionomial distributions». Computing 12 (3): 223–46.
Atkinson, AC, e MC Pearce. 1976. «The computer generation of beta, gamma and normal random variables». Journal of the Royal Statistical Society: Series A (General) 139 (4): 431–48.
Cheng, RCH. 1977. «The generation of gamma variables with non-integral shape parameter». Journal of the Royal Statistical Society: Series C (Applied Statistics) 26 (1): 71–75.
Härdle, Wolfgang Karl, Ostap Okhrin, e Yarema Okhrin. 2017. Basic elements of computational statistics. Springer.
Marsaglia, George. 1968. «Random numbers fall mainly in the planes». Proceedings of the National Academy of Sciences of the United States of America 61 (1): 25.
Neave, Henry R. 1973. «On using the Box-M üller transformation with multiplicative congruential pseudo-random number generators». Journal of the Royal Statistical Society: Series C (Applied Statistics) 22 (1): 92–97.
Wald, Abraham, e Jacob Wolfowitz. 1940. «On a test whether two samples are from the same population». The Annals of Mathematical Statistics 11 (2): 147–62.